Directory and doc rules

knitr::opts_chunk$set(
  echo = TRUE,         # Display code chunks
  eval = TRUE,         # Evaluate code chunks
  warning = FALSE,     # Hide warnings
  message = FALSE,     # Hide messages
  fig.width = 12,       # Set plot width in inches
  fig.height = 8,      # Set plot height in inches
  fig.align = "center" # Align plots to the center
)

Load packages

#install.packages(c("ggplot2", "sf", "viridis", "rnaturalearth", "rnaturalearthdata"))
library(tidyr)
#library(tidyverse)
library(ggplot2)
library(vegan)
library(sf)
library(viridis)
library(rnaturalearth)
library(rnaturalearthdata)

Load data

getwd()
## [1] "/Users/cmantegna/Documents/WDFWmussels/code"
#data has all sites, coordinates, p450, sod, condition factor, economic factor data
data<- read.csv("/Users/cmantegna/Documents/WDFWmussels/data/biomarkerfull.csv")
as.numeric("weight_initial", "weight_change", "length", "width", "height", "condition_factor", "economic_index")
## [1] NA

0’s and negative number management

SOD values that are negative are changed to 0 signifying that the value was undetectable.
p450 values that are 0 are missing, not actual 0’s.

#get rid of the p450 values that are 0
data <- data[data$p450 != 0, ]

#change any negative numbers to 0 to indicate not detectable
data$SOD <- ifelse(data$SOD <= 0, 0, data$SOD)

Map of Washington State

world <- ne_states(country = "united states of america", returnclass = "sf")
washington_map <- world[world$name == "Washington", ]
1
## [1] 1

p450 map

pmap<- ggplot() + 
  geom_sf(data = washington_map, fill = "lightgrey", color = "white") +
  geom_point(data = data, aes(x = longitude, y = latitude, color = p450), size = 1) +
  scale_color_viridis(option = "C", name = "p450") +
  theme_minimal() +
  labs(title = "Washington State - p450 Data")

print(pmap)

ggsave(plot=pmap, filename="/Users/cmantegna/Documents/WDFWmussels/output/p450map.png", width=15, height=8)

SOD map

smap<- ggplot() + 
  geom_sf(data = washington_map, fill = "lightgrey", color = "white") +
  geom_point(data = data, aes(x = longitude, y = latitude, color = SOD), size = 1) +
  scale_color_viridis(option = "C", name = "SOD") +
  theme_minimal() +
  labs(title = "Washington State - SOD Data")

print(smap)

ggsave(plot=smap, filename="/Users/cmantegna/Documents/WDFWmussels/output/sodmap.png", width=15, height=8)

Mapping both biomarkers

# Assuming washington_map is your sf object for Washington State
ggplot(data = washington_map) + 
  geom_sf(fill = "lightgrey", color = "white") +
  geom_point(data = data, aes(x = longitude, y = latitude, color = p450, size = SOD), alpha = 0.6) +
  scale_color_viridis(name = "p450", option = "D") +
  scale_size_continuous(name = "SOD", range = c(2, 6)) +
  theme_minimal() +
  labs(title = "Washington State: p450 and SOD")

Mapping with shapes to differenciate the biomarkers - change dataframe to accomodate

long_data <- pivot_longer(data, cols = c(SOD, p450), names_to = "type", values_to = "value")
ggplot(data = washington_map) + 
  geom_sf(fill = "lightgrey", color = "white") +
  geom_point(data = long_data, aes(x = longitude, y = latitude, color = value, shape = type), alpha = 0.6, size = 2) +
  scale_color_viridis(name = "value", option = "C") +
  scale_shape_manual(values = c("SOD" = 16, "p450" = 17)) + # 16: dot, 17: triangle
  theme_minimal() +
  labs(title = "Washington State: SOD & p450")

Mapping with Leaflet instead of ggplot; this is interactive and not helpful

if(!require(leaflet)) install.packages("leaflet")
library(leaflet)

leaflet(data) %>% addTiles() %>%
  addCircles(lng = ~longitude, lat = ~latitude, color = ~colorBin(c("red", "yellow", "green"), p450, bins = 3), radius = 500) %>%
  addCircles(lng = ~longitude, lat = ~latitude, color = ~colorBin(c("blue", "white", "black"), SOD, bins = 3), radius = 200)

Mapping with Plotly, also interactive and plotted on a grid that I don’t understand

if(!require(plotly)) install.packages("plotly")
library(plotly)

# Assuming data is your dataframe with lat, lon, p450, and SOD columns
fig <- plot_ly() %>%
  add_trace(data = data, x = ~longitude, y = ~latitude, type = 'scatter', mode = 'markers',
            marker = list(size = 10, color = ~p450, colorscale = 'Viridis', showscale = TRUE)) %>%
  add_trace(data = data, x = ~longitude, y = ~latitude, type = 'scatter', mode = 'markers',
            marker = list(size = 10, color = ~SOD, colorscale = 'Cividis', showscale = TRUE))

fig